Projektaufgaben Block 2

1. Nichtparametrisches Testen

Zwillingsstudie

Um zu testen, ob der Kindergartenbesuch einen signifikanten Einfluss auf die sozialen Fähigkeiten eines Kindes hat, führen wir einen zweiseitigen/einseitigen \(t\)-Test und einen zweiseitigen/einseitigen Wilcoxon-Vorzeichen-Test, jeweils zum Signifikanzniveau \(\alpha = 0.05\), durch.

# Enter data
x <- c(82,69,73,43,58,56,76,65)
y <- c(63,42,74,37,51,43,80,62)
testType <- c('two-sided', 'one-sided')
pValueT <- c(0,0); pValueWilcox <- c(0,0)

# t-test
pValueT[1] <- t.test(x, y, alternative = "two.sided", mu = 0, conf.level = 0.95, paired = TRUE)$p.value
pValueT[2] <- t.test(x, y, alternative = "greater", mu = 0, conf.level = 0.95, paired = TRUE)$p.value

# Wilcoxon signed rank test
pValueWilcox[1] <- wilcox.test(x, y, alternative = "two.sided", mu = 0, conf.level = 0.95, 
            paired = TRUE, conf.int = TRUE)$p.value
pValueWilcox[2] <- wilcox.test(x, y, alternative = "greater", mu = 0, conf.level = 0.95, 
            paired = TRUE, conf.int = TRUE)$p.value


knitr::kable(data.frame(testType, pValueT, pValueWilcox), caption = 'p-Werte')

p-Werte

testType	pValueT	pValueWilcox
two-sided	0.0489476	0.0546875
one-sided	0.0244738	0.0273438

Zweiseitige Tests

Wir können sehen, dass der \(t\)-Test die Nullhypothese ablehnt (\(p = 0.04895 < \alpha\)) und somit einen signifikanten Einfluss feststellt. Der Wilcoxon-Vorzeichen-Test hingegen (\(p = 0.0546875 > \alpha\)) lehnt die Nullhypothese nicht zum gegebenen Signifikanzniveau nicht ab.

Einseitige Tests

Im Fall zweiseitger Tests unter der gleichen Nullhypothese, und der Alternative dass die sozialen Fähigkeiten signifikant höher sind, lehnen nun beide Tests die Nullhypothese ab. Der Grund liegt darin, dass hier die \(p\)-Werte nur eine Ablehnregion in die Richtung “größer” abdecken müssen und aufgrund der Symmetrie der Tests deshalb nur halb so groß sind wie bei den zweiseitigen Tests. Beide Tests indizieren nun einen signifikanten Anstieg der sozialen Fähigkeiten durch den Kindergartenbesuch.

Testannahmen

Schließlich lohnt es, einen Blick auf die Testannahmen zu richten. Während der Wilcoxon-Test lediglich eine symmetrische Verteilung voraussetzt, nimmt der t-Test konkret die Normalverteilung, \(X_i - Y_i \sim N(0,\sigma^2)\), der Daten an. Da diese stärkere Annahme ihm auch zu seiner größeren Power verhilft, sollte man ihre Validität überprüfen, hier geschieht dies mit Hilfe eines QQ-Plot: Wie besonders am rechten Tail zu erkennen, kann die Normalverteilungsannahme durchaus in Frage gestellt werden.

Doch auch die Annahme symmetrischer Daten ist angesichts der “Schiefe” (2.25) der Daten anzuzweifeln, aber sie ist deutlich schwächer als die zusätzliche Annahme einer (bis auf Paramter) exakten Verteilung.

Das Hauptproblem dieser Stichprobe ist jedoch ihre geringe Größe von N = 8, jedwede Art von Statistik besitzt nur sehr geringe Aussagekraft.

\(t\)-Test vs. Wilcoxon-Vorzeichen-Test

Als nächstes schätzen wir die Power beider Tests in verschiedenen Settings ab. Dafür wird zunächst der wahre Wert $$ mit Fehlern verrauscht, die der Normal-/Cauchy-/ und Gleichverteilung folgen.

Blablabla

fnTestPowerMC <- function(fnError, n = 30, alpha = 0.05, nSim = 10^4, ...) {
  # This function estimates the probability of rejecting the null hypothesis of a
  # t-test and a Wilcoxon signed rank test using Monte Carlo simulations of iid
  # random variables X_i = \theta + \epsilon_i.
  # 
  # Args:
  #   fnError:  Function which generates random samples from a symmetric error 
  #             distribution \epsilon
  #   n:        Number of random samples used in tests
  #   alpha:    Significance level used in tests
  #   nSim:     Number of MC simulations of size n
  #   ...:      Further arguments to be passed to fnError
  #   
  # Returns:
  #   A list containing the following elements:
  #     $TProb:       MC estimation of rejection probability for the t-test
  #     $WilcoxProb:  MC estimation of rejection probability for the Wilcoxon
  #                   signed rank test
  
  # Perform MC simulation
  matX <- matrix(fnError(n*nSim, ...), ncol = n)
  
  # Define sub-functions which only return p-values from the two tests
  fnPvalT <- function(x) {
    return(t.test(x)$p.value)
  }
  fnPvalWilcox <- function(x) {
    return(wilcox.test(x)$p.value)
  }
  
  # Perform nSim number of tests with sample size n for each of the two tests
  vecPvalT <- apply(matX, 1, fnPvalT)
  vecPvalWilcox <- apply(matX, 1, fnPvalWilcox)
  
  # Compute and return estimations of rejection probabilities
  result <- list()
  result$TProb <- mean(vecPvalT < alpha)
  result$WilcoxProb <- mean(vecPvalWilcox < alpha)
  return(result)
}

# Set seed
set.seed(432)

# Normal errors
normalMean = c(0, 0.5, 1)
normalSD = c(1, 1, 1)
normalPowerT <- list()
normalPowerW <- list()

for (i in 1:3) {
  res <- fnTestPowerMC(fnError = rnorm, nSim = 10^4,
                       mean = normalMean[i], sd=normalSD[i])
  normalPowerT[i] <- res[1]
  normalPowerW[i] <- res[2]
}
# Cauchy errors (t-distribution with df = 1)
cauchyPowerT <- list()
cauchyPowerW <- list()
cauchyLoc = c(0, 0.5, 1)
cauchyScale = c(1,1,1)
for (i in 1:3) {
  res <- fnTestPowerMC(fnError = rcauchy, nSim = 10^4,
                       location = cauchyLoc[i], scale=cauchyScale[i])
  cauchyPowerT[i] <- res[1]
  cauchyPowerW[i] <- res[2]
}

# Uniform errors
uniMin <- c(-1, -0.5, 0)
uniMax <- c(1, 1.5, 2)
uniPowerT <- list()
uniPowerW <- list()

for (i in 1:3) {
  res <- fnTestPowerMC(fnError = runif, nSim = 10^4,
                       min = uniMin[i], max = uniMax[i])
  uniPowerT[i] <- res[1]
  uniPowerW[i] <- res[2]
}

normalverteiltes Rasuchen

mean	sd	power T-test	power Wilcox-test
0.0	1	0.0492	0.0506
0.5	1	0.7528	0.7387
1.0	1	0.9992	0.9990

uniformes Rauschen

mean	+- interval	power T-test	power Wilcox-test
0.0	1	0.0504	0.0495
0.5	1	0.9973	0.9871
1.0	1	1.0000	1.0000

cauchyverteiltes Rauschen

location	scale	power T-test	power Wilcox-test
0.0	1	0.0209	0.0527
0.5	1	0.0712	0.2928
1.0	1	0.2022	0.6927

Normalverteiltes Rauschen

Wie zu erwarten gibt es beim normalverteilten Rauschen kaum merkliche Unterschiede.

Uniformes Rauschen

t-test Überraschend gut im Vergleich zu Normal, quasi identisch. Liegt vermutlich an der

Cauchy-Rauschen

Deutlich schlechter als Wilcox, Cauchy-Verteilung stärker konzentriert als Normalverteilung, t-Test überschätzt Tails.

2. Dichteschätzung

Kerndichteschätzer

fnKernelDensityEst <- function(x, X, K, h) {
  # This function computes the kernel density estimates for a given sample X and
  # kernel K with bandwidth h at points x.
  # 
  # Args:
  #   x: Points for which the kernel density estimates are computed
  #   X: Data sample on which the kernel density estimation is fitted
  #   K: Kernel function to be used for smoothing
  #   h: Smoothing bandwidth
  #   
  # Returns:
  #   A vector of the kernel density estimates at points x
  
  # Compute kernel density estimation values using outer
  n <- length(X)
  mKernels <- K((1/h) * outer(X, x, FUN = "-"))
  return((1/(n*h)) * colSums(mKernels))
}

fnKernelPlot <- function(x, X, K, vh, title = NULL, trueDensity = NULL, ...) {
  # This function plots the kernel density estimates for different bandwidths.
  # If trueDensity is specified, the true density will be plotted as well.
  # 
  # Args:
  #   x:            Points for which the kernel density estimates are computed
  #   X:            Data sample on which the kernel density estimation is fitted
  #   K:            Kernel function to be used for smoothing
  #   vh:           Vector of smoothing bandwidths
  #   title:        Main title of the plot
  #   trueDensity:  True density to be plotted (e.g. dnorm)
  #   ...:          Further arguments passed to or from other methods
  #   
  # Returns:
  #   -
  
  nh <- length(vh)
  
  # Create color palette for number of bandwidths
  cols <- rainbow(nh, alpha = 1)
  
  # Initialize plot with first kernel density
  plot(x, fnKernelDensityEst(x, X, K, vh[1]), type = "l", col = cols[1], 
       main = title, xlab = "x", ylab = "Density", ...)
  
  # Add kernel density for each additional bandwidth
  if (nh > 1) {
    for (i in 2:nh) {
      lines(x, fnKernelDensityEst(x, X, K, vh[i]), col = cols[i])
    }
  }
  
  # Create legend
  legend <- sapply(vh, function(h) {h <- paste("h = ", round(h, 4))})
  
  # Plot true density if specified
  if (!is.null(trueDensity)) {
    lines(x, trueDensity(x), col = "black")
    legend <- c("True density", legend)
    cols <- c("black", cols)
  }
  
  # Plot legend
  legend("topright", legend = legend, col = cols, lty = 1, cex = 0.75)
}

# Define different smoothing kernels

fnRectangularKernel <- function(x) {
  return(0.5 * (abs(x) <= 1))
}

fnGaussianKernel <- function(x) {
  return((1/sqrt(2*pi)) * exp((-0.5) * x^2))
}

fnEpanechnikovKernel <- function(x) {
  return(0.75 * (1 - x^2) * (abs(x) <= 1))
}

Anwendung auf standardnormalverteilte Zufallszahlen

# Generate sample
set.seed(42)
n <- 50
X <- rnorm(n)

# Define points to estimate kernel density on
x <- seq(-4, 4, 0.01)

# Set different bandwidths
h <- c(0.1, bw.ucv(X), 1, 4)

## Warning in bw.ucv(X): minimum occurred at one end of the range

# Rectangular kernel density estimation
fnKernelPlot(x, X, fnRectangularKernel, h, 
             title = "Rectangular Kernel Density Estimation",
             trueDensity = dnorm, ylim = c(0, 0.6))

# Gaussian kernel density estimation
fnKernelPlot(x, X, fnGaussianKernel, h, 
             title = "Gaussian Kernel Density Estimation",
             trueDensity = dnorm, ylim = c(0, 0.6))

# Epanechnikov kernel density estimation
fnKernelPlot(x, X, fnGaussianKernel, h, 
             title = "Epanechnikov Kernel Density Estimation",
             trueDensity = dnorm, ylim = c(0, 0.6))

Anwendung auf gleichverteilte Zufallszahlen

# Generate sample
set.seed(42)
n <- 50
X <- runif(n)

# Define points to estimate kernel density on
x <- seq(0, 1, 0.01)

# Set different bandwidths
h <- c(bw.ucv(X), 0.1, 0.5)

# Rectangular kernel density estimation
fnKernelPlot(x, X, fnRectangularKernel, h, 
             title = "Rectangular Kernel Density Estimation",
             trueDensity = dunif, ylim = c(0, 2))

# Gaussian kernel density estimation
fnKernelPlot(x, X, fnGaussianKernel, h, 
             title = "Gaussian Kernel Density Estimation",
             trueDensity = dunif, ylim = c(0, 2))

# Epanechnikov kernel density estimation
fnKernelPlot(x, X, fnGaussianKernel, h, 
             title = "Epanechnikov Kernel Density Estimation",
             trueDensity = dunif, ylim = c(0, 2))

Anwendung auf Cauchy-verteilte Zufallszahlen

# Generate sample
set.seed(42)
n <- 50
X <- rt(n, df = 1)

# Define points to estimate kernel density on
x <- seq(-4, 4, 0.01)

# Set different bandwidths
h <- c(0.5, 1, 5, bw.ucv(X))

## Warning in bw.ucv(X): minimum occurred at one end of the range

# Define Cauchy density function
dcauchy <- function(x) {y <- dt(x, df = 1)}

# Rectangular kernel density estimation
fnKernelPlot(x, X, fnRectangularKernel, h, 
             title = "Rectangular Kernel Density Estimation",
             trueDensity = dcauchy, ylim = c(0, 0.4))

# Gaussian kernel density estimation
fnKernelPlot(x, X, fnGaussianKernel, h, 
             title = "Gaussian Kernel Density Estimation",
             trueDensity = dcauchy, ylim = c(0, 0.4))

# Epanechnikov kernel density estimation
fnKernelPlot(x, X, fnGaussianKernel, h, 
             title = "Epanechnikov Kernel Density Estimation",
             trueDensity = dcauchy, ylim = c(0, 0.4))

Anwendung auf faithful Datensatz

# Faithful kernel density estimation
X <- faithful$eruptions

# Define points to estimate kernel density on
x <- seq(min(X), max(X), 0.01)

# Set different bandwidths
h <- c(0.05, bw.ucv(X), 0.5, 1)

# Gaussian kernel density estimation
title <- "Gaussian Kernel Density Estimation for faithful eruption data"
fnKernelPlot(x, X, fnGaussianKernel, h, title = title)

Kreuzvalidierung zur Bandweitenwahl

fnJ <- function(X, K, h) {
  # This function computes the value of the unbiased cross validation criteria.
  # The minimum of this function in h gives an estimated “optimal” (in the
  # sense of mean integrated squared error) bandwidth.
  # 
  # Args:
  #   X: Data sample on which the kernel density estimation is fitted
  #   K: Kernel function to be used for smoothing
  #   h: Smoothing bandwidth
  #   
  # Returns:
  #   Function value
  
  n <- length(X)
  
  # Compute integral of squared kernel density estimator
  integrand <- function(x) {return(fnKernelDensityEst(x, X, K, h)^2)}
  SKDEInt <- integrate(integrand, -Inf, Inf)$value
  
  # Compute G
  mKernels <- K((1/h) * outer(X, X, FUN = "-"))
  diag(mKernels) <- 0
  G <- (1/(n*(n-1)*h)) * sum(mKernels)
  
  # Return function value of J
  return(SKDEInt - 2*G)
}

fnUCV <- function(X, K, hmin, hmax, tol = 0.1 * hmin) {
  # This function returns the minimum of the unbiased cross validation criteria
  # in h which is “optimal” in the sense of mean integrated squared error.
  # We use R's optimize function to find the minimum.
  # 
  # Args:
  #   X:    Data sample on which the kernel density estimation is fitted
  #   K:    Kernel function to be used for smoothing
  #   hmin: Lower bound for optimal h
  #   hmax: Upper bound for optimal h
  #   tol:  The desired accuracy
  #   
  # Returns:
  #   Optimal bandwidth h*
  
  # Find and return the minimum using optimize
  fnTarget <- function(h) {return(fnJ(X, K, h))}
  hOpt <- optimize(fnTarget, c(hmin, hmax), tol = tol)
  return(hOpt$minimum)
}

hmin <- 0.01
hmax <- 10
fnUCV(faithful$eruptions, fnGaussianKernel, hmin, hmax)

## [1] 0.102756

bw.ucv(faithful$eruptions)

## [1] 0.1019193

3. Bildentrauschen

# Load package
library(EBImage)

# Load image from parent directory
imgLenaColor <- readImage("../lena.png")

# Change image to grayscale
imgLenaGray <- channel(imgLenaColor, "gray")

# Display images
par(mfrow = c(1,2))
display(imgLenaColor, method = "raster")
display(imgLenaGray, method = "raster")

fnAddNoise <- function(img, rnoise, ...) {
  # This function adds noise specified by rnoise to a grayscale image. If the
  # image is not grayscale, it gets converted to grayscale first.
  # 
  # Args:
  #   img:    Image
  #   rnoise: Random noise generation function (e.g. rnorm for Gaussian noise)
  #   ...:    Further arguments passed to or from other methods
  #   
  # Returns:
  #   Grayscale image with added noise
  
  # Check if grayscale and convert if necessary
  if(colorMode(img) != 0) {img <- channel(img, "gray")}
  
  # Add noise
  m <- dim(img)[1]
  p <- dim(img)[2]
  imgNoise <- Image(imageData(img) + matrix(rnoise(m*p, ...), m, p))
  
  # Adjust values below 0 and above 1
  imgNoise[imgNoise < 0] <- 0
  imgNoise[imgNoise > 1] <- 1
  
  return(imgNoise)
}

imgLenaNoise1 <- fnAddNoise(imgLenaGray, rnorm, sd = 0.1)
imgLenaNoise2 <- fnAddNoise(imgLenaGray, rnorm, sd = 0.25)
imgLenaNoise3 <- fnAddNoise(imgLenaGray, rnorm, sd = 0.5)

par(mfrow = c(1,3))
display(imgLenaNoise1, method = "raster")
display(imgLenaNoise2, method = "raster")
display(imgLenaNoise3, method = "raster")

fnNadarayaWatson <- function(img, K, h) {
  # The Nadaraya-Watson-estimator with weights defined by kernel K and bandwidth
  # h for denoising/smoothing an image
  # 
  # Args:
  #   img:  Image to be denoised
  #   K:    Smoothing kernel used in weights
  #   h:    Bandwidth used in weights
  #   
  # Returns:
  #   Denoised/smoothed image
  
  # Check if grayscale and convert if necessary
  if(colorMode(img) != 0) {img <- channel(img, "gray")}
  
  # Get Y and dimensions
  Y <- imageData(img)
  m <- dim(img)[1]
  p <- dim(img)[2]
  
  # Generate M and N
  M <- K((1/h) * outer(1:m, 1:m, FUN = "-"))
  M <- M / rowSums(M)
  
  N <- K((1/h) * outer(1:p, 1:p, FUN = "-"))
  N <- N / rowSums(N)
  
  return(Image(M %*% Y %*% t(N)))
}

fnEvalNWGaussianNoise <- function(img, vsigma, vh) {
  # This function is a wrapper to compare the denoising results of images with 
  # different levels of Gaussian noise for the Nadaraya-Watson-estimator with
  # weights defined by the Gaussian kernel and rectangular kernel respectively 
  # for different bandwidths h.
  # 
  # Args:
  #   img:    Image
  #   sigma:  Vector of standard deviations used to add noise
  #   h:      Vector of bandwidths used in weights of the
  #           Nadaraya-Watson-estimator
  #   
  # Returns:
  #   -
  
  par(mfrow = c(length(vsigma), 1+2*length(vh)))
  
  for (sd in vsigma) {
    
    # Add noise to image and plot
    imgNoise <- fnAddNoise(img, rnorm, sd = sd)
    display(imgNoise, method = "raster")
    
    # NW-Denoising with Gaussian kernel
    for (h in vh) {
      display(fnNadarayaWatson(imgNoise, fnGaussianKernel, h), 
              method = "raster")
    }
    
    # NW-Denoising with rectangular kernel
    for (h in vh) {
      display(fnNadarayaWatson(imgNoise, fnRectangularKernel, h), 
              method = "raster")
    }
  }
}

sigma <- c(0.1, 0.25, 0.5)
h <- c(0.1, 1, 10)

# Lena
fnEvalNWGaussianNoise(imgLenaGray, sigma, h)

# Porsche
imgPorsche <- readImage("../porsche.jpeg")
fnEvalNWGaussianNoise(imgPorsche, sigma, h)

# 7. The weighted median is more robust than the mean